/*
 *  MWBounds.cpp
 *  T3nsors2
 *
 *  Created by Michael Barriault on 10-07-02.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "MW.h"
#include <cmath>
#include <iostream>
using namespace t3;
using std::cout;
using std::endl;

void MW::Bounds(scalar* Phi, scalar* Psi, vector* E, vector* B, bool xp, bool yp, bool zp) {
	double RT2 = sqrt(2.);
	double RT3 = sqrt(3.);
	
	scalar* phi = scalar::init(O,"Phi");
	scalar* psi = scalar::init(O,"Psi");
	vector* e = vector::init(O,3,"E");
	vector* b = vector::init(O,3,"B");

	FOR(o,this->O->size()) {
		phi->at(o) = Phi->at(o);
		psi->at(o) = Psi->at(o);
		FOR(a,3) e->c(a)->at(o) = E->c(a)->at(o);
		FOR(a,3) b->c(a)->at(o) = B->c(a)->at(o);
	}
	
	if ( not xp ) FOR(j,O->at(1)) FOR(k,O->at(2)) { // x=0,1 plane
		Phi->pt(0,j,k) = 0.5*(phi->pt(0,j,k)+e->c(0)->pt(0,j,k));
		Psi->pt(0,j,k) = 0.5*(psi->pt(0,j,k)+b->c(0)->pt(0,j,k));
		E->c(0)->pt(0,j,k) = 0.5*(e->c(0)->pt(0,j,k)+phi->pt(0,j,k));
		E->c(1)->pt(0,j,k) = 0.5*(e->c(1)->pt(0,j,k)-b->c(2)->pt(0,j,k));
		E->c(2)->pt(0,j,k) = 0.5*(e->c(2)->pt(0,j,k)+b->c(1)->pt(0,j,k));
		B->c(0)->pt(0,j,k) = 0.5*(b->c(0)->pt(0,j,k)+psi->pt(0,j,k));
		B->c(1)->pt(0,j,k) = 0.5*(b->c(1)->pt(0,j,k)+e->c(2)->pt(0,j,k));
		B->c(2)->pt(0,j,k) = 0.5*(b->c(2)->pt(0,j,k)-e->c(1)->pt(0,j,k));
		Phi->pt(-1,j,k) = 0.5*(phi->pt(-1,j,k)-e->c(0)->pt(-1,j,k));
		Psi->pt(-1,j,k) = 0.5*(psi->pt(-1,j,k)-b->c(0)->pt(-1,j,k));
		E->c(0)->pt(-1,j,k) = 0.5*(e->c(0)->pt(-1,j,k)-phi->pt(-1,j,k));
		E->c(1)->pt(-1,j,k) = 0.5*(e->c(1)->pt(-1,j,k)+b->c(2)->pt(-1,j,k));
		E->c(2)->pt(-1,j,k) = 0.5*(e->c(2)->pt(-1,j,k)-b->c(1)->pt(-1,j,k));
		B->c(0)->pt(-1,j,k) = 0.5*(b->c(0)->pt(-1,j,k)-psi->pt(-1,j,k));
		B->c(1)->pt(-1,j,k) = 0.5*(b->c(1)->pt(-1,j,k)-e->c(2)->pt(-1,j,k));
		B->c(2)->pt(-1,j,k) = 0.5*(b->c(2)->pt(-1,j,k)+e->c(1)->pt(-1,j,k));
	}
	
	if ( not yp ) FOR(i,O->at(0)) FOR(k,O->at(2)) { // y=0,1 plane
		Phi->pt(i,0,k) = 0.5*(phi->pt(i,0,k)+e->c(1)->pt(i,0,k));
		Psi->pt(i,0,k) = 0.5*(psi->pt(i,0,k)+b->c(1)->pt(i,0,k));
		E->c(0)->pt(i,0,k) = 0.5*(e->c(0)->pt(i,0,k)+b->c(2)->pt(i,0,k));
		E->c(1)->pt(i,0,k) = 0.5*(e->c(1)->pt(i,0,k)+phi->pt(i,0,k));
		E->c(2)->pt(i,0,k) = 0.5*(e->c(2)->pt(i,0,k)-b->c(0)->pt(i,0,k));
		B->c(0)->pt(i,0,k) = 0.5*(b->c(0)->pt(i,0,k)-e->c(2)->pt(i,0,k));
		B->c(1)->pt(i,0,k) = 0.5*(b->c(1)->pt(i,0,k)+psi->pt(i,0,k));
		B->c(2)->pt(i,0,k) = 0.5*(b->c(2)->pt(i,0,k)+e->c(0)->pt(i,0,k));
		Phi->pt(i,-1,k) = 0.5*(phi->pt(i,-1,k)-e->c(1)->pt(i,-1,k));
		Psi->pt(i,-1,k) = 0.5*(psi->pt(i,-1,k)-b->c(1)->pt(i,-1,k));
		E->c(0)->pt(i,-1,k) = 0.5*(e->c(0)->pt(i,-1,k)-b->c(2)->pt(i,-1,k));
		E->c(1)->pt(i,-1,k) = 0.5*(e->c(1)->pt(i,-1,k)-phi->pt(i,-1,k));
		E->c(2)->pt(i,-1,k) = 0.5*(e->c(2)->pt(i,-1,k)+b->c(0)->pt(i,-1,k));
		B->c(0)->pt(i,-1,k) = 0.5*(b->c(0)->pt(i,-1,k)+e->c(2)->pt(i,-1,k));
		B->c(1)->pt(i,-1,k) = 0.5*(b->c(1)->pt(i,-1,k)-psi->pt(i,-1,k));
		B->c(2)->pt(i,-1,k) = 0.5*(b->c(2)->pt(i,-1,k)-e->c(0)->pt(i,-1,k));
	}
	
	if ( not zp ) FOR(i,O->at(0)) FOR(j,O->at(1)) { // z=0,1 plane
		Phi->pt(i,j,0) = 0.5*(phi->pt(i,j,0)+e->c(2)->pt(i,j,0));
		Psi->pt(i,j,0) = 0.5*(psi->pt(i,j,0)+b->c(2)->pt(i,j,0));
		E->c(0)->pt(i,j,0) = 0.5*(e->c(0)->pt(i,j,0)-b->c(1)->pt(i,j,0));
		E->c(1)->pt(i,j,0) = 0.5*(e->c(1)->pt(i,j,0)+b->c(0)->pt(i,j,0));
		E->c(2)->pt(i,j,0) = 0.5*(e->c(2)->pt(i,j,0)+phi->pt(i,j,0));
		b->c(0)->pt(i,j,0) = 0.5*(b->c(0)->pt(i,j,0)+e->c(1)->pt(i,j,0));
		b->c(1)->pt(i,j,0) = 0.5*(b->c(1)->pt(i,j,0)-e->c(0)->pt(i,j,0));
		b->c(2)->pt(i,j,0) = 0.5*(b->c(2)->pt(i,j,0)+psi->pt(i,j,0));
		Phi->pt(i,j,-1) = 0.5*(phi->pt(i,j,-1)-e->c(2)->pt(i,j,-1));
		Psi->pt(i,j,-1) = 0.5*(psi->pt(i,j,-1)-b->c(2)->pt(i,j,-1));
		E->c(0)->pt(i,j,-1) = 0.5*(e->c(0)->pt(i,j,-1)+b->c(1)->pt(i,j,-1));
		E->c(1)->pt(i,j,-1) = 0.5*(e->c(1)->pt(i,j,-1)-b->c(0)->pt(i,j,-1));
		E->c(2)->pt(i,j,-1) = 0.5*(e->c(2)->pt(i,j,-1)-phi->pt(i,j,-1));
		B->c(0)->pt(i,j,-1) = 0.5*(b->c(0)->pt(i,j,-1)-e->c(1)->pt(i,j,-1));
		B->c(1)->pt(i,j,-1) = 0.5*(b->c(1)->pt(i,j,-1)+e->c(0)->pt(i,j,-1));
		B->c(2)->pt(i,j,-1) = 0.5*(b->c(2)->pt(i,j,-1)-psi->pt(i,j,-1));
	}
	
	
	if ( not xp and not yp ) FOR(k,O->at(2)) { // x,y = constant lines
		Phi->pt(0,0,k) = 0.5*(phi->pt(0,0,k)+(e->c(0)->pt(0,0,k)+e->c(1)->pt(0,0,k))/RT2);
		Psi->pt(0,0,k) = 0.5*(psi->pt(0,0,k)+(b->c(0)->pt(0,0,k)+b->c(1)->pt(0,0,k))/RT2);
		E->c(0)->pt(0,0,k) = 0.5*(e->c(0)->pt(0,0,k)+(phi->pt(0,0,k)+b->c(2)->pt(0,0,k))/RT2);
		E->c(1)->pt(0,0,k) = 0.5*(e->c(1)->pt(0,0,k)+(phi->pt(0,0,k)-b->c(2)->pt(0,0,k))/RT2);
		E->c(2)->pt(0,0,k) = 0.5*(e->c(2)->pt(0,0,k)+(b->c(1)->pt(0,0,k)-b->c(0)->pt(0,0,k))/RT2);
		B->c(0)->pt(0,0,k) = 0.5*(b->c(0)->pt(0,0,k)+(psi->pt(0,0,k)-e->c(2)->pt(0,0,k))/RT2);
		B->c(1)->pt(0,0,k) = 0.5*(b->c(1)->pt(0,0,k)+(psi->pt(0,0,k)+e->c(2)->pt(0,0,k))/RT2);
		B->c(2)->pt(0,0,k) = 0.5*(b->c(2)->pt(0,0,k)+(e->c(0)->pt(0,0,k)-e->c(1)->pt(0,0,k))/RT2);
		Phi->pt(-1,-1,k) = 0.5*(phi->pt(-1,-1,k)-(e->c(0)->pt(-1,-1,k)+e->c(1)->pt(-1,-1,k))/RT2);
		Psi->pt(-1,-1,k) = 0.5*(psi->pt(-1,-1,k)-(b->c(0)->pt(-1,-1,k)+b->c(1)->pt(-1,-1,k))/RT2);
		E->c(0)->pt(-1,-1,k) = 0.5*(e->c(0)->pt(-1,-1,k)-(phi->pt(-1,-1,k)+b->c(2)->pt(-1,-1,k))/RT2);
		E->c(1)->pt(-1,-1,k) = 0.5*(e->c(1)->pt(-1,-1,k)-(phi->pt(-1,-1,k)-b->c(2)->pt(-1,-1,k))/RT2);
		E->c(2)->pt(-1,-1,k) = 0.5*(e->c(2)->pt(-1,-1,k)+(b->c(0)->pt(-1,-1,k)-b->c(1)->pt(-1,-1,k))/RT2);
		B->c(0)->pt(-1,-1,k) = 0.5*(b->c(0)->pt(-1,-1,k)-(psi->pt(-1,-1,k)-e->c(2)->pt(-1,-1,k))/RT2);
		B->c(1)->pt(-1,-1,k) = 0.5*(b->c(1)->pt(-1,-1,k)-(psi->pt(-1,-1,k)+e->c(2)->pt(-1,-1,k))/RT2);
		B->c(2)->pt(-1,-1,k) = 0.5*(b->c(2)->pt(-1,-1,k)+(e->c(1)->pt(-1,-1,k)-e->c(0)->pt(-1,-1,k))/RT2);
		Phi->pt(-1,0,k) = 0.5*(phi->pt(-1,0,k)+(e->c(1)->pt(-1,0,k)-e->c(0)->pt(-1,0,k))/RT2);
		Psi->pt(-1,0,k) = 0.5*(psi->pt(-1,0,k)+(b->c(1)->pt(-1,0,k)-b->c(0)->pt(-1,0,k))/RT2);
		E->c(0)->pt(-1,0,k) = 0.5*(e->c(0)->pt(-1,0,k)-(phi->pt(-1,0,k)-b->c(2)->pt(-1,0,k))/RT2);
		E->c(1)->pt(-1,0,k) = 0.5*(e->c(1)->pt(-1,0,k)+(phi->pt(-1,0,k)+b->c(2)->pt(-1,0,k))/RT2);
		E->c(2)->pt(-1,0,k) = 0.5*(e->c(2)->pt(-1,0,k)-(b->c(0)->pt(-1,0,k)+b->c(1)->pt(-1,0,k))/RT2);
		B->c(0)->pt(-1,0,k) = 0.5*(b->c(0)->pt(-1,0,k)-(psi->pt(-1,0,k)+e->c(2)->pt(-1,0,k))/RT2);
		B->c(1)->pt(-1,0,k) = 0.5*(b->c(1)->pt(-1,0,k)+(psi->pt(-1,0,k)-e->c(2)->pt(-1,0,k))/RT2);
		B->c(2)->pt(-1,0,k) = 0.5*(b->c(2)->pt(-1,0,k)+(e->c(0)->pt(-1,0,k)+e->c(1)->pt(-1,0,k))/RT2);
		Phi->pt(0,-1,k) = 0.5*(phi->pt(0,-1,k)+(e->c(0)->pt(0,-1,k)-e->c(1)->pt(0,-1,k))/RT2);
		Psi->pt(0,-1,k) = 0.5*(psi->pt(0,-1,k)+(b->c(0)->pt(0,-1,k)-b->c(1)->pt(0,-1,k))/RT2);
		E->c(0)->pt(0,-1,k) = 0.5*(e->c(0)->pt(0,-1,k)+(phi->pt(0,-1,k)-b->c(2)->pt(0,-1,k))/RT2);
		E->c(1)->pt(0,-1,k) = 0.5*(e->c(1)->pt(0,-1,k)-(phi->pt(0,-1,k)+b->c(2)->pt(0,-1,k))/RT2);
		E->c(2)->pt(0,-1,k) = 0.5*(e->c(2)->pt(0,-1,k)+(b->c(0)->pt(0,-1,k)+b->c(1)->pt(0,-1,k))/RT2);
		B->c(0)->pt(0,-1,k) = 0.5*(b->c(0)->pt(0,-1,k)+(psi->pt(0,-1,k)+e->c(2)->pt(0,-1,k))/RT2);
		B->c(1)->pt(0,-1,k) = 0.5*(b->c(1)->pt(0,-1,k)-(psi->pt(0,-1,k)-e->c(2)->pt(0,-1,k))/RT2);
		B->c(2)->pt(0,-1,k) = 0.5*(b->c(2)->pt(0,-1,k)-(e->c(0)->pt(0,-1,k)+e->c(1)->pt(0,-1,k))/RT2);
	}

	if ( not yp and not zp ) FOR(i,O->at(0)) { // y,z = constant lines
		Phi->pt(i,0,0) = 0.5*(phi->pt(i,0,0)+(e->c(1)->pt(i,0,0)+e->c(2)->pt(i,0,0))/RT2);
		Psi->pt(i,0,0) = 0.5*(psi->pt(i,0,0)+(b->c(1)->pt(i,0,0)+b->c(2)->pt(i,0,0))/RT2);
		E->c(0)->pt(i,0,0) = 0.5*(e->c(0)->pt(i,0,0)+(b->c(2)->pt(i,0,0)-b->c(1)->pt(i,0,0))/RT2);
		E->c(1)->pt(i,0,0) = 0.5*(e->c(1)->pt(i,0,0)+(phi->pt(i,0,0)+b->c(0)->pt(i,0,0))/RT2);
		E->c(2)->pt(i,0,0) = 0.5*(e->c(2)->pt(i,0,0)+(phi->pt(i,0,0)-b->c(0)->pt(i,0,0))/RT2);
		B->c(0)->pt(i,0,0) = 0.5*(b->c(0)->pt(i,0,0)+(e->c(1)->pt(i,0,0)-e->c(2)->pt(i,0,0))/RT2);
		B->c(1)->pt(i,0,0) = 0.5*(b->c(1)->pt(i,0,0)+(psi->pt(i,0,0)-e->c(0)->pt(i,0,0))/RT2);
		B->c(2)->pt(i,0,0) = 0.5*(b->c(2)->pt(i,0,0)+(psi->pt(i,0,0)+e->c(0)->pt(i,0,0))/RT2);
		Phi->pt(i,-1,-1) = 0.5*(phi->pt(i,-1,-1)-(e->c(1)->pt(i,-1,-1)+e->c(2)->pt(i,-1,-1))/RT2);
		Psi->pt(i,-1,-1) = 0.5*(psi->pt(i,-1,-1)-(b->c(1)->pt(i,-1,-1)+b->c(2)->pt(i,-1,-1))/RT2);
		E->c(0)->pt(i,-1,-1) = 0.5*(e->c(0)->pt(i,-1,-1)+(b->c(1)->pt(i,-1,-1)-b->c(2)->pt(i,-1,-1))/RT2);
		E->c(1)->pt(i,-1,-1) = 0.5*(e->c(1)->pt(i,-1,-1)-(phi->pt(i,-1,-1)+b->c(0)->pt(i,-1,-1))/RT2);
		E->c(2)->pt(i,-1,-1) = 0.5*(e->c(2)->pt(i,-1,-1)-(phi->pt(i,-1,-1)-b->c(0)->pt(i,-1,-1))/RT2);
		B->c(0)->pt(i,-1,-1) = 0.5*(b->c(0)->pt(i,-1,-1)+(e->c(2)->pt(i,-1,-1)-e->c(1)->pt(i,-1,-1))/RT2);
		B->c(1)->pt(i,-1,-1) = 0.5*(b->c(1)->pt(i,-1,-1)-(psi->pt(i,-1,-1)-e->c(0)->pt(i,-1,-1))/RT2);
		B->c(2)->pt(i,-1,-1) = 0.5*(b->c(2)->pt(i,-1,-1)-(psi->pt(i,-1,-1)+e->c(0)->pt(i,-1,-1))/RT2);
		Phi->pt(i,-1,0) = 0.5*(phi->pt(i,-1,0)+(e->c(2)->pt(i,-1,0)-e->c(1)->pt(i,-1,0))/RT2);
		Psi->pt(i,-1,0) = 0.5*(psi->pt(i,-1,0)+(b->c(2)->pt(i,-1,0)-b->c(1)->pt(i,-1,0))/RT2);
		E->c(0)->pt(i,-1,0) = 0.5*(e->c(0)->pt(i,-1,0)-(b->c(1)->pt(i,-1,0)+b->c(2)->pt(i,-1,0))/RT2);
		E->c(1)->pt(i,-1,0) = 0.5*(e->c(1)->pt(i,-1,0)-(phi->pt(i,-1,0)-b->c(0)->pt(i,-1,0))/RT2);
		E->c(2)->pt(i,-1,0) = 0.5*(e->c(2)->pt(i,-1,0)+(phi->pt(i,-1,0)+b->c(0)->pt(i,-1,0))/RT2);
		B->c(0)->pt(i,-1,0) = 0.5*(b->c(0)->pt(i,-1,0)+(e->c(1)->pt(i,-1,0)+e->c(2)->pt(i,-1,0))/RT2);
		B->c(1)->pt(i,-1,0) = 0.5*(b->c(1)->pt(i,-1,0)-(psi->pt(i,-1,0)+e->c(0)->pt(i,-1,0))/RT2);
		B->c(2)->pt(i,-1,0) = 0.5*(b->c(2)->pt(i,-1,0)+(psi->pt(i,-1,0)-e->c(0)->pt(i,-1,0))/RT2);
		Phi->pt(i,0,-1) = 0.5*(phi->pt(i,0,-1)+(e->c(1)->pt(i,0,-1)-e->c(2)->pt(i,0,-1))/RT2);
		Psi->pt(i,0,-1) = 0.5*(psi->pt(i,0,-1)+(b->c(1)->pt(i,0,-1)-b->c(2)->pt(i,0,-1))/RT2);
		E->c(0)->pt(i,0,-1) = 0.5*(e->c(0)->pt(i,0,-1)+(b->c(1)->pt(i,0,-1)+b->c(2)->pt(i,0,-1))/RT2);
		E->c(1)->pt(i,0,-1) = 0.5*(e->c(1)->pt(i,0,-1)+(phi->pt(i,0,-1)-b->c(0)->pt(i,0,-1))/RT2);
		E->c(2)->pt(i,0,-1) = 0.5*(e->c(2)->pt(i,0,-1)-(phi->pt(i,0,-1)+b->c(0)->pt(i,0,-1))/RT2);
		B->c(0)->pt(i,0,-1) = 0.5*(b->c(0)->pt(i,0,-1)-(e->c(1)->pt(i,0,-1)+e->c(2)->pt(i,0,-1))/RT2);
		B->c(1)->pt(i,0,-1) = 0.5*(b->c(1)->pt(i,0,-1)+(psi->pt(i,0,-1)+e->c(0)->pt(i,0,-1))/RT2);
		B->c(2)->pt(i,0,-1) = 0.5*(b->c(2)->pt(i,0,-1)-(psi->pt(i,0,-1)-e->c(0)->pt(i,0,-1))/RT2);
	}
	
	if ( not xp and not zp ) FOR(j,O->at(1)) { // x,z = constant lines
		Phi->pt(0,j,0) = 0.5*(phi->pt(0,j,0)+(e->c(2)->pt(0,j,0)+e->c(0)->pt(0,j,0))/RT2);
		Psi->pt(0,j,0) = 0.5*(psi->pt(0,j,0)+(b->c(2)->pt(0,j,0)+b->c(0)->pt(0,j,0))/RT2);
		E->c(0)->pt(0,j,0) = 0.5*(e->c(0)->pt(0,j,0)+(phi->pt(0,j,0)-b->c(1)->pt(0,j,0))/RT2);
		E->c(1)->pt(0,j,0) = 0.5*(e->c(1)->pt(0,j,0)+(b->c(0)->pt(0,j,0)-b->c(2)->pt(0,j,0))/RT2);
		E->c(2)->pt(0,j,0) = 0.5*(e->c(2)->pt(0,j,0)+(phi->pt(0,j,0)+b->c(1)->pt(0,j,0))/RT2);
		B->c(0)->pt(0,j,0) = 0.5*(b->c(0)->pt(0,j,0)+(psi->pt(0,j,0)+e->c(1)->pt(0,j,0))/RT2);
		B->c(1)->pt(0,j,0) = 0.5*(b->c(1)->pt(0,j,0)+(e->c(2)->pt(0,j,0)-e->c(0)->pt(0,j,0))/RT2);
		B->c(2)->pt(0,j,0) = 0.5*(b->c(2)->pt(0,j,0)+(psi->pt(0,j,0)-e->c(1)->pt(0,j,0))/RT2);
		Phi->pt(-1,j,-1) = 0.5*(phi->pt(-1,j,-1)-(e->c(2)->pt(-1,j,-1)+e->c(0)->pt(-1,j,-1))/RT2);
		Psi->pt(-1,j,-1) = 0.5*(psi->pt(-1,j,-1)-(b->c(2)->pt(-1,j,-1)+b->c(0)->pt(-1,j,-1))/RT2);
		E->c(0)->pt(-1,j,-1) = 0.5*(e->c(0)->pt(-1,j,-1)-(phi->pt(-1,j,-1)-b->c(1)->pt(-1,j,-1))/RT2);
		E->c(1)->pt(-1,j,-1) = 0.5*(e->c(1)->pt(-1,j,-1)+(b->c(2)->pt(-1,j,-1)-b->c(0)->pt(-1,j,-1))/RT2);
		E->c(2)->pt(-1,j,-1) = 0.5*(e->c(2)->pt(-1,j,-1)-(phi->pt(-1,j,-1)+b->c(1)->pt(-1,j,-1))/RT2);
		B->c(0)->pt(-1,j,-1) = 0.5*(b->c(0)->pt(-1,j,-1)-(psi->pt(-1,j,-1)+e->c(1)->pt(-1,j,-1))/RT2);
		B->c(1)->pt(-1,j,-1) = 0.5*(b->c(1)->pt(-1,j,-1)+(e->c(0)->pt(-1,j,-1)-e->c(2)->pt(-1,j,-1))/RT2);
		B->c(2)->pt(-1,j,-1) = 0.5*(b->c(2)->pt(-1,j,-1)-(psi->pt(-1,j,-1)-e->c(1)->pt(-1,j,-1))/RT2);
		Phi->pt(0,j,-1) = 0.5*(phi->pt(0,j,-1)+(e->c(0)->pt(0,j,-1)-e->c(2)->pt(0,j,-1))/RT2);
		Psi->pt(0,j,-1) = 0.5*(psi->pt(0,j,-1)+(b->c(0)->pt(0,j,-1)-b->c(2)->pt(0,j,-1))/RT2);
		E->c(0)->pt(0,j,-1) = 0.5*(e->c(0)->pt(0,j,-1)+(phi->pt(0,j,-1)+b->c(1)->pt(0,j,-1))/RT2);
		E->c(1)->pt(0,j,-1) = 0.5*(e->c(1)->pt(0,j,-1)-(b->c(2)->pt(0,j,-1)+b->c(0)->pt(0,j,-1))/RT2);
		E->c(2)->pt(0,j,-1) = 0.5*(e->c(2)->pt(0,j,-1)-(phi->pt(0,j,-1)-b->c(1)->pt(0,j,-1))/RT2);
		B->c(0)->pt(0,j,-1) = 0.5*(b->c(0)->pt(0,j,-1)+(psi->pt(0,j,-1)-e->c(1)->pt(0,j,-1))/RT2);
		B->c(1)->pt(0,j,-1) = 0.5*(b->c(1)->pt(0,j,-1)+(e->c(2)->pt(0,j,-1)+e->c(0)->pt(0,j,-1))/RT2);
		B->c(2)->pt(0,j,-1) = 0.5*(b->c(2)->pt(0,j,-1)-(psi->pt(0,j,-1)+e->c(1)->pt(0,j,-1))/RT2);
		Phi->pt(-1,j,0) = 0.5*(phi->pt(-1,j,0)+(e->c(2)->pt(-1,j,0)-e->c(0)->pt(-1,j,0))/RT2);
		Psi->pt(-1,j,0) = 0.5*(psi->pt(-1,j,0)+(b->c(2)->pt(-1,j,0)-b->c(0)->pt(-1,j,0))/RT2);
		E->c(0)->pt(-1,j,0) = 0.5*(e->c(0)->pt(-1,j,0)-(phi->pt(-1,j,0)+b->c(1)->pt(-1,j,0))/RT2);
		E->c(1)->pt(-1,j,0) = 0.5*(e->c(1)->pt(-1,j,0)+(b->c(2)->pt(-1,j,0)+b->c(0)->pt(-1,j,0))/RT2);
		E->c(2)->pt(-1,j,0) = 0.5*(e->c(2)->pt(-1,j,0)+(phi->pt(-1,j,0)-b->c(1)->pt(-1,j,0))/RT2);
		B->c(0)->pt(-1,j,0) = 0.5*(b->c(0)->pt(-1,j,0)-(psi->pt(-1,j,0)-e->c(1)->pt(-1,j,0))/RT2);
		B->c(1)->pt(-1,j,0) = 0.5*(b->c(1)->pt(-1,j,0)-(e->c(2)->pt(-1,j,0)+e->c(0)->pt(-1,j,0))/RT2);
		B->c(2)->pt(-1,j,0) = 0.5*(b->c(2)->pt(-1,j,0)+(psi->pt(-1,j,0)+e->c(1)->pt(-1,j,0))/RT2);
	}
	
	if ( not xp and not yp and not zp ) {
		{ // point 0,0,0
			Phi->pt(0,0,0) = 0.5*(phi->pt(0,0,0)+(e->c(0)->pt(0,0,0)+e->c(1)->pt(0,0,0)+e->c(2)->pt(0,0,0))/RT3);
			Psi->pt(0,0,0) = 0.5*(psi->pt(0,0,0)+(b->c(0)->pt(0,0,0)+b->c(1)->pt(0,0,0)+b->c(2)->pt(0,0,0))/RT3);
			E->c(0)->pt(0,0,0) = 0.5*(e->c(0)->pt(0,0,0)+(phi->pt(0,0,0)-b->c(1)->pt(0,0,0)+b->c(2)->pt(0,0,0))/RT3);
			E->c(1)->pt(0,0,0) = 0.5*(e->c(1)->pt(0,0,0)+(phi->pt(0,0,0)+b->c(0)->pt(0,0,0)-b->c(2)->pt(0,0,0))/RT3);
			E->c(2)->pt(0,0,0) = 0.5*(e->c(2)->pt(0,0,0)+(phi->pt(0,0,0)-b->c(0)->pt(0,0,0)+b->c(1)->pt(0,0,0))/RT3);
			B->c(0)->pt(0,0,0) = 0.5*(b->c(0)->pt(0,0,0)+(psi->pt(0,0,0)+e->c(1)->pt(0,0,0)-e->c(2)->pt(0,0,0))/RT3);
			B->c(1)->pt(0,0,0) = 0.5*(b->c(1)->pt(0,0,0)+(psi->pt(0,0,0)-e->c(0)->pt(0,0,0)+e->c(2)->pt(0,0,0))/RT3);
			B->c(2)->pt(0,0,0) = 0.5*(b->c(2)->pt(0,0,0)+(psi->pt(0,0,0)+e->c(0)->pt(0,0,0)-e->c(1)->pt(0,0,0))/RT3);
		}
		{ // point 0,0,1
			Phi->pt(0,0,-1) = 0.5*(phi->pt(0,0,-1)+(e->c(0)->pt(0,0,-1)+e->c(1)->pt(0,0,-1)-e->c(2)->pt(0,0,-1))/RT3);
			Psi->pt(0,0,-1) = 0.5*(psi->pt(0,0,-1)+(b->c(0)->pt(0,0,-1)+b->c(1)->pt(0,0,-1)-b->c(2)->pt(0,0,-1))/RT3);
			E->c(0)->pt(0,0,-1) = 0.5*(e->c(0)->pt(0,0,-1)+(phi->pt(0,0,-1)+b->c(1)->pt(0,0,-1)+b->c(2)->pt(0,0,-1))/RT3);
			E->c(1)->pt(0,0,-1) = 0.5*(e->c(1)->pt(0,0,-1)+(phi->pt(0,0,-1)-b->c(0)->pt(0,0,-1)-b->c(2)->pt(0,0,-1))/RT3);
			E->c(2)->pt(0,0,-1) = 0.5*(e->c(2)->pt(0,0,-1)-(phi->pt(0,0,-1)+b->c(0)->pt(0,0,-1)-b->c(1)->pt(0,0,-1))/RT3);
			B->c(0)->pt(0,0,-1) = 0.5*(b->c(0)->pt(0,0,-1)+(psi->pt(0,0,-1)-e->c(1)->pt(0,0,-1)-e->c(2)->pt(0,0,-1))/RT3);
			B->c(1)->pt(0,0,-1) = 0.5*(b->c(1)->pt(0,0,-1)+(psi->pt(0,0,-1)+e->c(0)->pt(0,0,-1)+e->c(2)->pt(0,0,-1))/RT3);
			B->c(2)->pt(0,0,-1) = 0.5*(b->c(2)->pt(0,0,-1)-(psi->pt(0,0,-1)-e->c(0)->pt(0,0,-1)+e->c(1)->pt(0,0,-1))/RT3);
		}
		{ // point 0,1,0
			Phi->pt(0,-1,0) = 0.5*(phi->pt(0,-1,0)+(e->c(0)->pt(0,-1,0)-e->c(1)->pt(0,-1,0)+e->c(2)->pt(0,-1,0))/RT3);
			Psi->pt(0,-1,0) = 0.5*(psi->pt(0,-1,0)+(b->c(0)->pt(0,-1,0)-b->c(1)->pt(0,-1,0)+b->c(2)->pt(0,-1,0))/RT3);
			E->c(0)->pt(0,-1,0) = 0.5*(e->c(0)->pt(0,-1,0)+(phi->pt(0,-1,0)-b->c(1)->pt(0,-1,0)-b->c(2)->pt(0,-1,0))/RT3);
			E->c(1)->pt(0,-1,0) = 0.5*(e->c(1)->pt(0,-1,0)-(phi->pt(0,-1,0)-b->c(0)->pt(0,-1,0)+b->c(2)->pt(0,-1,0))/RT3);
			E->c(2)->pt(0,-1,0) = 0.5*(e->c(2)->pt(0,-1,0)+(phi->pt(0,-1,0)+b->c(0)->pt(0,-1,0)+b->c(1)->pt(0,-1,0))/RT3);
			B->c(0)->pt(0,-1,0) = 0.5*(b->c(0)->pt(0,-1,0)+(psi->pt(0,-1,0)+e->c(1)->pt(0,-1,0)+e->c(2)->pt(0,-1,0))/RT3);
			B->c(1)->pt(0,-1,0) = 0.5*(b->c(1)->pt(0,-1,0)-(psi->pt(0,-1,0)+e->c(0)->pt(0,-1,0)-e->c(2)->pt(0,-1,0))/RT3);
			B->c(2)->pt(0,-1,0) = 0.5*(b->c(2)->pt(0,-1,0)+(psi->pt(0,-1,0)-e->c(0)->pt(0,-1,0)-e->c(1)->pt(0,-1,0))/RT3);
		}
		{ // point 1,0,0
			Phi->pt(-1,0,0) = 0.5*(phi->pt(-1,0,0)+(-e->c(0)->pt(-1,0,0)+e->c(1)->pt(-1,0,0)+e->c(2)->pt(-1,0,0))/RT3);
			Psi->pt(-1,0,0) = 0.5*(psi->pt(-1,0,0)+(-b->c(0)->pt(-1,0,0)+b->c(1)->pt(-1,0,0)+b->c(2)->pt(-1,0,0))/RT3);
			E->c(0)->pt(-1,0,0) = 0.5*(e->c(0)->pt(-1,0,0)-(phi->pt(-1,0,0)+b->c(1)->pt(-1,0,0)-b->c(2)->pt(-1,0,0))/RT3);
			E->c(1)->pt(-1,0,0) = 0.5*(e->c(1)->pt(-1,0,0)+(phi->pt(-1,0,0)+b->c(0)->pt(-1,0,0)+b->c(2)->pt(-1,0,0))/RT3);
			E->c(2)->pt(-1,0,0) = 0.5*(e->c(2)->pt(-1,0,0)+(phi->pt(-1,0,0)-b->c(0)->pt(-1,0,0)-b->c(1)->pt(-1,0,0))/RT3);
			B->c(0)->pt(-1,0,0) = 0.5*(b->c(0)->pt(-1,0,0)-(psi->pt(-1,0,0)-e->c(1)->pt(-1,0,0)+e->c(2)->pt(-1,0,0))/RT3);
			B->c(1)->pt(-1,0,0) = 0.5*(b->c(1)->pt(-1,0,0)+(psi->pt(-1,0,0)-e->c(0)->pt(-1,0,0)-e->c(2)->pt(-1,0,0))/RT3);
			B->c(2)->pt(-1,0,0) = 0.5*(b->c(2)->pt(-1,0,0)+(psi->pt(-1,0,0)+e->c(0)->pt(-1,0,0)+e->c(1)->pt(-1,0,0))/RT3);
		}
		{ // point 1,1,0
			Phi->pt(-1,-1,0) = 0.5*(phi->pt(-1,-1,0)+(-e->c(0)->pt(-1,-1,0)-e->c(1)->pt(-1,-1,0)+e->c(2)->pt(-1,-1,0))/RT3);
			Psi->pt(-1,-1,0) = 0.5*(psi->pt(-1,-1,0)+(-b->c(0)->pt(-1,-1,0)-b->c(1)->pt(-1,-1,0)+b->c(2)->pt(-1,-1,0))/RT3);
			E->c(0)->pt(-1,-1,0) = 0.5*(e->c(0)->pt(-1,-1,0)-(phi->pt(-1,-1,0)+b->c(1)->pt(-1,-1,0)+b->c(2)->pt(-1,-1,0))/RT3);
			E->c(1)->pt(-1,-1,0) = 0.5*(e->c(1)->pt(-1,-1,0)-(phi->pt(-1,-1,0)-b->c(0)->pt(-1,-1,0)-b->c(2)->pt(-1,-1,0))/RT3);
			E->c(2)->pt(-1,-1,0) = 0.5*(e->c(2)->pt(-1,-1,0)+(phi->pt(-1,-1,0)+b->c(0)->pt(-1,-1,0)-b->c(1)->pt(-1,-1,0))/RT3);
			B->c(0)->pt(-1,-1,0) = 0.5*(b->c(0)->pt(-1,-1,0)-(psi->pt(-1,-1,0)-e->c(1)->pt(-1,-1,0)-e->c(2)->pt(-1,-1,0))/RT3);
			B->c(1)->pt(-1,-1,0) = 0.5*(b->c(1)->pt(-1,-1,0)-(psi->pt(-1,-1,0)+e->c(0)->pt(-1,-1,0)+e->c(2)->pt(-1,-1,0))/RT3);
			B->c(2)->pt(-1,-1,0) = 0.5*(b->c(2)->pt(-1,-1,0)+(psi->pt(-1,-1,0)-e->c(0)->pt(-1,-1,0)+e->c(1)->pt(-1,-1,0))/RT3);
		}
		{ // point 1,0,1
			Phi->pt(-1,0,-1) = 0.5*(phi->pt(-1,0,-1)+(-e->c(0)->pt(-1,0,-1)+e->c(1)->pt(-1,0,-1)-e->c(2)->pt(-1,0,-1))/RT3);
			Psi->pt(-1,0,-1) = 0.5*(psi->pt(-1,0,-1)+(-b->c(0)->pt(-1,0,-1)+b->c(1)->pt(-1,0,-1)-b->c(2)->pt(-1,0,-1))/RT3);
			E->c(0)->pt(-1,0,-1) = 0.5*(e->c(0)->pt(-1,0,-1)-(phi->pt(-1,0,-1)-b->c(1)->pt(-1,0,-1)-b->c(2)->pt(-1,0,-1))/RT3);
			E->c(1)->pt(-1,0,-1) = 0.5*(e->c(1)->pt(-1,0,-1)+(phi->pt(-1,0,-1)-b->c(0)->pt(-1,0,-1)+b->c(2)->pt(-1,0,-1))/RT3);
			E->c(2)->pt(-1,0,-1) = 0.5*(e->c(2)->pt(-1,0,-1)-(phi->pt(-1,0,-1)+b->c(0)->pt(-1,0,-1)+b->c(1)->pt(-1,0,-1))/RT3);
			B->c(0)->pt(-1,0,-1) = 0.5*(b->c(0)->pt(-1,0,-1)-(psi->pt(-1,0,-1)+e->c(1)->pt(-1,0,-1)+e->c(2)->pt(-1,0,-1))/RT3);
			B->c(1)->pt(-1,0,-1) = 0.5*(b->c(1)->pt(-1,0,-1)+(psi->pt(-1,0,-1)+e->c(0)->pt(-1,0,-1)-e->c(2)->pt(-1,0,-1))/RT3);
			B->c(2)->pt(-1,0,-1) = 0.5*(b->c(2)->pt(-1,0,-1)-(psi->pt(-1,0,-1)-e->c(0)->pt(-1,0,-1)-e->c(1)->pt(-1,0,-1))/RT3);
		}
		{ // point 0,1,1
			Phi->pt(0,-1,-1) = 0.5*(phi->pt(0,-1,-1)+(e->c(0)->pt(0,-1,-1)-e->c(1)->pt(0,-1,-1)-e->c(2)->pt(0,-1,-1))/RT3);
			Psi->pt(0,-1,-1) = 0.5*(psi->pt(0,-1,-1)+(b->c(0)->pt(0,-1,-1)-b->c(1)->pt(0,-1,-1)-b->c(2)->pt(0,-1,-1))/RT3);
			E->c(0)->pt(0,-1,-1) = 0.5*(e->c(0)->pt(0,-1,-1)+(phi->pt(0,-1,-1)+b->c(1)->pt(0,-1,-1)-b->c(2)->pt(0,-1,-1))/RT3);
			E->c(1)->pt(0,-1,-1) = 0.5*(e->c(1)->pt(0,-1,-1)-(phi->pt(0,-1,-1)+b->c(0)->pt(0,-1,-1)+b->c(2)->pt(0,-1,-1))/RT3);
			E->c(2)->pt(0,-1,-1) = 0.5*(e->c(2)->pt(0,-1,-1)-(phi->pt(0,-1,-1)-b->c(0)->pt(0,-1,-1)-b->c(1)->pt(0,-1,-1))/RT3);
			B->c(0)->pt(0,-1,-1) = 0.5*(b->c(0)->pt(0,-1,-1)+(psi->pt(0,-1,-1)-e->c(1)->pt(0,-1,-1)+e->c(2)->pt(0,-1,-1))/RT3);
			B->c(1)->pt(0,-1,-1) = 0.5*(b->c(1)->pt(0,-1,-1)-(psi->pt(0,-1,-1)-e->c(0)->pt(0,-1,-1)-e->c(2)->pt(0,-1,-1))/RT3);
			B->c(2)->pt(0,-1,-1) = 0.5*(b->c(2)->pt(0,-1,-1)-(psi->pt(0,-1,-1)+e->c(0)->pt(0,-1,-1)+e->c(1)->pt(0,-1,-1))/RT3);
		}
		{ // point 1,1,1
			Phi->pt(-1,-1,-1) = 0.5*(phi->pt(-1,-1,-1)-(e->c(0)->pt(-1,-1,-1)+e->c(1)->pt(-1,-1,-1)+e->c(2)->pt(-1,-1,-1))/RT3);
			Psi->pt(-1,-1,-1) = 0.5*(psi->pt(-1,-1,-1)-(b->c(0)->pt(-1,-1,-1)+b->c(1)->pt(-1,-1,-1)+b->c(2)->pt(-1,-1,-1))/RT3);
			E->c(0)->pt(-1,-1,-1) = 0.5*(e->c(0)->pt(-1,-1,-1)-(phi->pt(-1,-1,-1)-b->c(1)->pt(-1,-1,-1)+b->c(2)->pt(-1,-1,-1))/RT3);
			E->c(1)->pt(-1,-1,-1) = 0.5*(e->c(1)->pt(-1,-1,-1)-(phi->pt(-1,-1,-1)+b->c(0)->pt(-1,-1,-1)-b->c(2)->pt(-1,-1,-1))/RT3);
			E->c(2)->pt(-1,-1,-1) = 0.5*(e->c(2)->pt(-1,-1,-1)-(phi->pt(-1,-1,-1)-b->c(0)->pt(-1,-1,-1)+b->c(1)->pt(-1,-1,-1))/RT3);
			B->c(0)->pt(-1,-1,-1) = 0.5*(b->c(0)->pt(-1,-1,-1)-(psi->pt(-1,-1,-1)+e->c(1)->pt(-1,-1,-1)-e->c(2)->pt(-1,-1,-1))/RT3);
			B->c(1)->pt(-1,-1,-1) = 0.5*(b->c(1)->pt(-1,-1,-1)-(psi->pt(-1,-1,-1)-e->c(0)->pt(-1,-1,-1)+e->c(2)->pt(-1,-1,-1))/RT3);
			B->c(2)->pt(-1,-1,-1) = 0.5*(b->c(2)->pt(-1,-1,-1)-(psi->pt(-1,-1,-1)+e->c(0)->pt(-1,-1,-1)-e->c(1)->pt(-1,-1,-1))/RT3);
		}
	}
	
	if ( xp ) FOR(j,O->at(1)) FOR(k,O->at(2)) {
		Phi->pt(0,j,k) = Phi->pt(-1,j,k) = 0.5*(phi->pt(0,j,k)+phi->pt(-1,j,k));
		Psi->pt(0,j,k) = Psi->pt(-1,j,k) = 0.5*(psi->pt(0,j,k)+psi->pt(-1,j,k));
		FOR(a,3) E->c(a)->pt(0,j,k) = E->c(a)->pt(-1,j,k) = 0.5*(e->c(a)->pt(0,j,k)+e->c(a)->pt(-1,j,k));
		FOR(a,3) B->c(a)->pt(0,j,k) = B->c(a)->pt(-1,j,k) = 0.5*(b->c(a)->pt(0,j,k)+b->c(a)->pt(-1,j,k));
	}
	if ( yp ) FOR(i,O->at(0)) FOR(k,O->at(2)) {
		Phi->pt(i,0,k) = Phi->pt(i,-1,k) = 0.5*(phi->pt(i,0,k)+phi->pt(i,-1,k));
		Psi->pt(i,0,k) = Psi->pt(i,-1,k) = 0.5*(psi->pt(i,0,k)+psi->pt(i,-1,k));
		FOR(a,3) E->c(a)->pt(i,0,k) = E->c(a)->pt(i,-1,k) = 0.5*(e->c(a)->pt(i,0,k)+e->c(a)->pt(i,-1,k));
		FOR(a,3) B->c(a)->pt(i,0,k) = B->c(a)->pt(i,-1,k) = 0.5*(b->c(a)->pt(i,0,k)+b->c(a)->pt(i,-1,k));
	}
	if ( zp ) FOR(i,O->at(0)) FOR(j,O->at(1)) {
		Phi->pt(i,j,0) = Phi->pt(i,j,-1) = 0.5*(phi->pt(i,j,0)+phi->pt(i,j,-1));
		Psi->pt(i,j,0) = Psi->pt(i,j,-1) = 0.5*(psi->pt(i,j,0)+psi->pt(i,j,-1));
		FOR(a,3) E->c(a)->pt(i,j,0) = E->c(a)->pt(i,j,-1) = 0.5*(e->c(a)->pt(i,j,0)+e->c(a)->pt(i,j,-1));
		FOR(a,3) B->c(a)->pt(i,j,0) = B->c(a)->pt(i,j,-1) = 0.5*(b->c(a)->pt(i,j,0)+b->c(a)->pt(i,j,-1));
	}
	
	
/*	FOR(o,this->O->size()) {
		phi->at(o) = Phi->at(o);
		psi->at(o) = Psi->at(o);
		FOR(a,3) e->c(a)->at(o) = E->c(a)->at(o);
		FOR(a,3) b->c(a)->at(o) = B->c(a)->at(o);
	}*/
	
	phi->release(); psi->release(); e->release(); b->release();
	
}
